IOP* 
' E SCHOOL 
MONTEREY.  CALIFORNIA   aiHAQ 


NPS-63Wu76101 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


A  FINITE  ELEMENT  PREDICTION  MODEL 

WITH  VARIABLE  ELEMENT 

SIZES 

by 

R.  G.  Kelley,  Jr. 

R,  T.  Williams 

October  1976 

Techn 

ical 

Report  Period:   January  - 

September 

1976 

Aiproved  for  public  release;  distribution  unlimited. 

FEDDOCS  apared  for:   Naval  Environmental  Prediction  Research  Facility, 

D208.14/2:NPS-63Wu76101  Monterey,  California 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 

Rear  Admiral  Isham  W.  Linder  Jack  R.  Borsting 

Superintendent  Provost 

ABSTRACT 

There  are  a  variety  of  meteorological  forecast  problems  which  require 
high  spatial  resolution  in  only  a  limited  area.   An  important  example  of 
this  type  of  problem  is  the  prediction  of  tropical  cyclones.   This  study 
tests  a  simple  finite  element  prediction  model  with  a  variable  element 
size.   The  shallow  water  equations  are  used  and  the  motion  is  confined  in 
a  periodic  channel  on  a  f-plane.   The  Galerkin  technique  is  applied  to 
linear  basis  functions  on  triangular  elements.   The  model  uses  leapfrog 
time  differencing  and  periodic  restarts.   The  model  is  tested  with  a  wave 
imbedded  in  a  mean  flow  and  also  with  an  isolated  vortex.   The  experiments 
with  a  uniform  element  size  show  excellent  phase  propagation,  but  some 
small  scale  noise  is  generated.   The  introduction  of  momentum  diffusion 
terms  helps  to  control  the  noise.   The  model  is  also  tested  with  elements 
which  decrease  abruptly  in  scale  along  a  line  and  with  elements  which  de- 
crease smoothly.   Both  of  these  cases  generate  more  noise  than  with  uniform 
elements. 
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ABSTRACT 

There  are  a  variety  of  meteorological  forecast  problems  which  require 
high  spatial  resolution  in  only  a  limited  area.   An  important  example  of 
this  type  of  problem  is  the  prediction  of  tropical  cyclones.   This  study 
tests  a  simple  finite  element  prediction  model  with  a  variable  element 
size.   The  shallow  water  equations  are  used  and  the  motion  is  confined  in 
a  periodic  channel  on  a  f-plane.   The  Galerkin  technique  is  applied  to 
linear  basis  functions  on  triangular  elements.   The  model  uses  leapfrog 
time  differencing  and  periodic  restarts.   The  model  is  tested  with  a  wave 
imbedded  in  a  mean  flow  and  also  with  an  isolated  vortex.   The  experiments 
with  a  uniform  element  size  show  excellent  phase  propagation,  but  some 
small  scale  noise  is  generated.   The  introduction  of  momentum  diffusion 
terms  helps  to  control  the  noise.   The  model  is  also  tested  with  elements 
which  decrease  abruptly  in  scale  along  a  line  and  with  elements  which  de- 
crease smoothly.   Both  of  these  cases  generate  more  noise  than  with 
uniform  elements. 
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LIST  OF  SYMBOLS  AND  ABBREVIATIONS 

A  area 

a  radius  of  earth  or  an  elementary  coordinate 

a  zonal  wind  component 

3  meridional  wind  component 

b  elementary  coordinate  or  a  vector 

d  a  length  along  a  diagonal 

d  a  constant  length  along  a  diagonal 

e  element  index 

f  Coriolis  parameter 

fQ  Coriolis  parameter  at  central  latitude 

(J)  geopotential  height 

g  test  function 

y  geopotential  height 

i  nodal  point  index 

j  nodal  point  index 

K  an  arbitrary  diagonal  factor 

K^  diffusion  coefficient 

k  nodal  point  index 

a£  {    }    differential  operator 

L  area  coordinate 

LFM  limited  fine  mesh 

1  nodal  point  index 

A  longitude 

mb  millibar 

N  length  of  an  array 

n  a  time  level  or  a  Gauss-Seidel  iteration 

Q,  angular  velocity  of  earth 

r  radial  distance  from  vortex  center 

r  maximum  radial  extent  of  vortex 

t  time  or  a  test  function 

6  latitude 

0  central  latitude 
o 

U  mean  zonal  wind 

u  zonal  wind  component 

V  basis  or  test  function 

VT  tangential  wind 

v  meridional  wind  component 

W  channel  width 

x  longitudinal  coordinate 

Xt  channel  length 

y  latitudinal  coordinate 

z  a  field  variable 

<   >    inner  product  (area  integration) 

[    ]     matrix 

{    }    vector 

an  approximation 


I.   INTRODUCTION 

Operational  numerical  models,  because  of  computational  limitations, 
are  able  to  depict  and  predict  atmospheric  processes  only  down  to  a 
limited  scale.   Until  the  introduction  of  the  National  Meteorological 
Center's  limited  fine  mesh  (LFM)  model,  which  is  superimposed  on  the 
coarse  mesh  Northern  Hemisphere  stereographic  grid,  grids  were  uniform 
on  any  given  map  projection  and  the  advent  of  LFM  has  brought  increased 
accuracy  in  wave  phase  and  amplitude  prediction  due  to  smaller  trunca- 
tion errors  [Houghton  and  Irvine  (1976)].   Nevertheless,  problems  arise 
at  the  boundaries  of  the  fine  mesh  grid  because  of  the  abrupt  change  in 
grid  size  and  the  need  to  furnish  boundary  values  interpolated  from  the 
coarse  mesh.   Hence  it  seems  logical  that  a  scheme  which  permits  a 
gradual  change  in  resolution  would  be  desirable. 

A  relatively  new  technique,  known  as  the  finite  element  method  (FEM) , 
seemed  suited  to  the  problem  of  a  non-uniform  grid.   In  theory  the  de- 
pendent variable  may  be  predicted  by  this  method  at  any  number  of  loca- 
tions chosen  at  will.   The  main  thrust  of  this  particular  investigation 
was  therefore  to  test  this  hypothesis.   Secondary  emphasis  was  placed 
on  determining  the  computational  limitations,  computer  time  and  storage 
requirements  of  such  a  procedure. 

Only  recently  have  the  atmospheric  sciences  seen  the  use  of  the 
technique.   Among  the  notable  contributions  are  those  of  Cullen  (1973, 
1974,  1976);  Lee,  Gresho,  and  Sani  (personal  communication);  and  Hinsman 
(1975) .   Since  the  finite  element  method  (FEM)  is  an  implicit  scheme  in 


space,  physically  reasonable  solutions  involve  solving  a  large  system 
of  linear  equations.   This  implies  substantial  memory  requirements  and, 
depending  on  the  technique  of  matrix  inversion,  possibly  considerable 
CPU  time. 

Specifically,  the  technique  of  Galerkin  (1915)  was  used  to  transform 
the  equation  set  into  an  integrable  form.   Galerkin' s  method  has  recent- 
ly been  proven  invaluable  in  engineering  disciplines.   For  example,  the 
field  of  indeterminate  structural  analysis,  pressure  vessel  and  air 
frame  design,  and  continuum  mechanics  have  all  greatly  benefited. 

The  equation  set  used  was  the  "shallow  water",  barotropic  equations 
which  approximately  describe  atmospheric  motion.  These  equations  model  a 
simple,  physically  reasonable  form  of  atmospheric  motion.   Additional 
simplification  was  gained  by  assuming  the  coriolous  parameter,  f ,  was 
constant  and  equal  to  its  value  at  the  mid-channel  latitude.   This  assump- 
tion does  not  degrade  model  performance  significantly  since  the  actual 
value  of  f  varies  slowly  over  the  domain  studied.   The  value  of  this 
approximation  lies  principally  in  reducing  the  number  of  terms  in  the 
equations  by  two  and  eliminating  the  "Beta  effect".   With  no  f  variation, 
sinusoidal  waves  can  move  in  an  east-west  direction  without  any  resultant 
latitudinal  tilt. 

Two  basic  types  of  domains  and  initial  conditions  were  tested.  Both 
domains  were  square  with  a  channel  length  of  approximately  3500  kilometers. 
The  channel  width  was  equivalent  to  a  region  extending  from  the  equator 
to  30  degrees  north  latitude.  The  channel  was  periodic  in  the  east-west 
direction. 


Simple  experiments  were  carried  out  for  sinusoidal  initial  conditions 
which  consisted  of  a  single  wave  in  the  longitudinal  or  x-direction  and 
a  half  wave  in  the  latitudinal  or  y-direction.   A  geopotential  surface 
tilt  was  added  in  the  y-direction  to  allow  for  any  mean  zonal  flow  imposed. 
Numerous  tests  were  performed  with  these  initial  conditions  using  various 
domain  subdivisions.   All  domain  subdivisions  were  done  in  a  regular  manner 
by  varying  the  grid  spacing  in  either  or  both  the  x  and  y  directions  over 
the  domain.   Additional  tests  were  also  conducted  on  some  of  these  meshes 
with  an  analytic  vortex. 

Final  and  more  complex  tests  involved  a  domain  on  which  the  grid  reso- 
lution decreased  regularly  toward  the  center  of  the  domain.   Initial  condi- 
tions used  consisted  of  those  mentioned  above  and  a  localized  disturbance 
or  vortex  which  spanned  the  central  region  of  varying  resolution. 

The  simpler  tests  were  undertaken  to  determine  the  effects  of  various 
smooth  and  abrupt  grid  changes  on  model  performance.   The  final,  complex 
model  was  evaluated  with  an  eye  towards  ultimate  implementation  as  an 
operational  tool  for  forecasting  tropical  cyclone  movement.   It  was  hoped 
that  the  flexibility  of  choosing  a  varying  grid  resolution  would  be  a 
significant  improvement  over  current  finite  difference  models.   It  seems 
desirable  to  have  a  scheme  by  which  one  can  vary  grid  resolution  smoothly 
to  accommodate  the  scale  of  motion  changes  that  occur  proceeding  radially 
inward  into  a  typhoon. 
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II.   FINITE  LINEAR  ELEMENTS 

As  is  true  with  the  spectral  method,  the  FEM  assumes  that  the  depen- 
dent variable  can  be  represented  by  a  known  function  of  space.   The  two 
techniques  differ  since  spectral  approximations  are  continuous  and  are 
infinitely  dif ferentiable  while  FEM  functions  are  defined  very  locally 
(over  an  element) ,  vanishing  elsewhere,  with  limited  continuity.   Gener- 
ally, FEM  functions  are  low  order  polynomials,  defined  as  planar  surfaces 
locally  in  the  case  of  linear  elements.   That  is,  for  a  given  variable, 
$    ,  one  assumes  that  over  each  element  the  basis  or  local  support  function 

V.   has  the  form:   V.  =  a  +  bx  +  cy  where  a,  b,  and  c  are  constants. 

3  3 

Figure  1  shows  an  exploded,  three  dimensional  view  of  a  small  domain 
which  has  been  subdivided  into  triangular  elements.   Above  this  domain 
lies  the  group  of  linear  basis  functions  which  are  pyramids  rising  to  an 
approximation  of  (f)  over  a  nodal  point.   Above  these  pyramids  is  the  seg- 
mented, planar  approximation  of  (J>  with  dashed  cross-sections  of  the  exact 
<j)  surface  at  evenly  spaced  Ax  values.   Note  that  $     is  shown  varying 
only  in  one  direction  in  this  case  and  that  the  <p     approximations,   $  , 
are  shown  to  be  exact  at  the  nodal  points.   The  latter  fact  is,  of 
course,  not  true  in  general.   What  is  required  for  piecewise  continuity 
in  (J) .   is  that  no  nodal  point  may  be  located  on  the  side  of  a  triangle. 

Once  a  form  for  the  unknown  as  been  assumed,  Galerkin's  technique  may 
be  applied.   In  its  basic  form,  this  method  assumes  an  unknown  has  the 
following  approximate  form  locally  over  a  nodal  point,  j  and  its  surround- 
ing elements: 


(p.    =   Y.V.  (II-l) 

3  3    3 
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Figure  1.   Linear  finite  element  approximation. 
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(with  the  repeated  subscript  implying  a  sum  over  that  subscript) 

V.   is  defined  as  the  basis  or  support  function  and  has  maximum  value 

of  unity  at  the  node,  sloping  to  zero  at  surrounding  elements'  outer 

boundaries.   Thus  with  the  aid  of  Figure  1  it  is  seen  the  Y .  =  d> .   at 

3  J 

any  nodal  point  j . 

To  solve  a  differential  equation  in  ty    ,  assume  (J)   is  of  the  form 
(II-l) ,  multiply  by  a  suitable  "test"  function,  and  integrate  over  the 
domain.   These  test  functions,  if  chosen  wisely,  greatly  simplify  the 
resulting  equation  and  make  it  easily  integrable.   A  generalized  example 
treats,  in  one  dimension,  the  equation  (with  c&   {   }   a  differential 
operator) : 

<*f{(t>}  -  f(x)  =0      0  <  x  <  a 


An  approximate  form  for  (J)   is  assumed,   (J)  ;  the  test  function, 

g(x),  is  applied,  and  the  integration  performed.   The  resulting  equation 

is: 

a 
f      [^{$1  -  f(x)]g(x)  dx  =  0  (II-2) 

o 

Since   [  o£ {(J)}-f  (x)  ]  =  R  where  R  is  the  error  in  the  approximation, 
this  integration  should  minimize  R,  with  g(x)  properly  chosen.   Galerkin's 
choice  for  g  (x)  is  the  same  form  as  that  for  (J).  ,  (II-l),  with  the  co- 
efficient Y-  =  1-   The  attraction  of  this  approach  is  that  since  both 
g.   and  <J) .   are  defined  only  locally,  i.e.  over  elements  surrounding 
nodal  point  j,  the  result  is  a  linear  matrix  equation  which  is  banded. 
This  results  from  a  partial  integration  of  (II-2) .   Another  useful  feature 
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is  that  linear  elements  have  domain  integrations  which  are  analytically 
Obtainable. 

This  direct  integration  is  made  straightforward  through  the  use  of 
a  natural  or  area  coordinate  scheme  for  each  element  as  explained  in 
Desai  and  Abel  (1972) .   In  two  space  dimensions  the  location  of  a  point 
P(x,y)  inside  an  element  (see  Figure  2)  can  be  defined  by  describing  the 
three  areas:   A  ,  A  ,  and  A  ,  subtended  by  lines  from  P  to  the  triangle's 
vertices  (x . ,  y.,  j=l,  2,  3).   Zienkiewicz  (1971)  shows  how  the  location 
of  P  can  be  transformed  into  area  coordinates  via: 


x  =  L  .x  . 

]  : 

y  =  L.y.  (II-3) 

3    3 

1  =  L   +  L2  +  L3 


with:     j  =  1,  2,  3 


L.=  A. /A 
J         3 


A  =  total  elementary  area 
If  the  differential  equation  to  be  solved  is: 


||  ♦$  ♦♦!,„>  =o 


xn  <  x  <  x^  (II-4) 

1  —   —   2 


Galerkin's  technique  yields 

y2   x2 

Y-   /    /   (v^  ""  v^„  +  Vn.  )V_,dxdy  =  ( 


xx     ly     l   j 


Yl   Xl 


with  variable  subscript  implying  partial  differentiation. 
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Y 

i 


-*~X 


X3*3 


P^i.kJL3) 


2*2 


Figure   2.      Area   coordinates, 
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This  matrix  equation  may  be  solved  for  the  y .   coefficients  once 
the  area  integration  is  performed.   Zienkiewicz  show  that  for  triangular 
elements  with  linear  basis  and  test  functions  a  general  area  integration 
formula  is: 


// 


_  m  „  n  ,  ,     m!n! 
Ll  L2  dxdy  =  (m+n+2) i 


2A 


(II-5) 


It  is  noted  that,  for  m,  n>l,  non-linear  terms  are  being  formulated. 

For  example,  with  m=l  and  n=2,  the  value  of  (II-5)  becomes  —  . 

3  0 

Transforming  equation  set  (II-3)  into  matrix  form  and  solving  for 
the  L  vector  one  obtains: 


\ 


1 
2A 


2A 


2A 


b   a 
1   1 


2A    b2   a2 


b3   a3 


(ii-e: 


with  a. ,  b.  as  defined  in  Figure  3. 

Differentiation  of  (II-6)  shows  the  differential  operations  desired 
to  complete  the  integration  of  (4)  are: 


8 
dx  ' 

i 
2A 

9 

3L. 

i 

9 

a. 

i 

8 

3y  " 

2A 

9L. 

i 

(Henceforth  a  repeated  subscript  will  imply  a  sum  from  one  through 
three.) 
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Y 

A 


i 


b2=Vyi 


i 


a  *  x    -  x 
3      2        1 


T 


Vv 


-1- 

b3-yry2 


■**-x 


Figure  3.   Local  coordinate  system. 
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Applying  the  chain  rule  in  natural  coordinates; 


3v.    b,   3v.    b„   3v.    b„   9v. 
1    1    1    2    i    3    1 


3x   '  2A   3L     2A   3L2    2A   9L_ 


but 


3v. 

-5-i  =0  for  i  ^  k 

3Lk 

3v. 

and  since   V.  =  L,  ;  -kt —  =  1   for   i 
1     1     dJj, 
k 


so 


3V.    b.  3v.    a. 

li       ,      li 


and    -r —  -  — -  (II-7) 


3x    2A  3y    2A 


Therefore  equation  (II-4)  is  easily  integrated  on  an  elementary  level, 
These  results  can  then  be  distributed  throughout  the  resultant  coeffi- 
cient  matrix  and  the  solution  vector   {<(>}   obtained  via  a  linear  equation 
solver. 
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III.   EQUATIONS  AND  BOUNDARY  CONDITIONS 

A  form  of  the  barotropic,  shallow  water,  primitive  equations  was 
selected  which  was  developed  by  Phillips  (1959) .   The  original  intent 
of  this  study  was  to  develop  a  complete,  barotropic  model  for  a  tropical 
channel.   Time  limitations  forced  the  use  of  a  model  which  was  purely 
cartesian.   That  is,  although  the  original  equation  set  was  designed 
in  mercator  coordinates,  a  constant  map  factor  of  unity  was  used.   In 
theory,  and  based  on  other  work  done  by  this  author,  extension  of  the 
Galerkin  formulation  of  the  equation  set  to  include  complete  sphericity 
is  easily  done. 

The  set  in  cartesian  coordinates  becomes 


!£  =  -[-!-  w  +  t-  wn                                (iii-D 

dt  dx         dy 

-r—  =  -br~  +  u  -r—  +  v  -r— ]  +  2Qv   sin  6                      (III-2) 

dt  dx      dx      dy             o 

77-  =  -[t^-   +  u  tc—  +  v  7T- ]    -    2Slu  sin  9                      (III-3) 

dt  dy      dx      dy              o 


with  the  transform  from  spherical  to  the  cartesian  system  accomplished 
using: 

x  =  aA 

cos  9 

y  =  a  In 


1-sin 

and  9  =  latitude 

A  =  longitude 
9   -  15° 
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The  only  mathematical  boundary  conditions  imposed  were  those  of  no  cross 
channel  flow  at  the  latitudinal  boundaries  and  of  cyclic  continuity  in 
the  longitudinal  direction.   Through  the  course  of  model  development 
other  conditions  were  tried  at  the  boundaries  due  to  physical  constraints, 

A  geos trophic  balance  at  channel  walls  was  imposed  in  the  mass  continuity 

1 
equation ( I II-l ) ,  i.e.: 


|£  =  -  2ftu  sin  9  (III-4) 

dy  o 


Additionally  it  was  thought  necessary  to  impose  a  zero  normal 
gradient  of  zonal  wind  at  the  walls.   During  some  early  experiments 
with  an  abruptly  varying  element  size,  it  was  found  that  imposing  this 
condition  on  the  zonal  equation  of  motion  yielded  instability  after  six 
to  12  hours  of  time  integration.   Later  experimentation  showed  the  zonal 
equation  required  no  explicit  boundary  conditions. 


This  condition  is  consistent  with  that  of  no  cross  boundary  flow 
and  follows  from  inserting  V=0  in  equation  (III-3) . 
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IV.   DIFFUSION 

Early  experiments  with  a  coarse  mesh  indicated  that  some  type  of 

filtering  would  be  required  since  short  waves  of  twice  the  coarse  grid 

- 
size  were  readily  generated.   These  short  waves  were  evident  in  the  dis- 
turbance patterns  of  the  geopotential  and  zonal  wind  fields.   Since  ex- 
perience and  knowledge  of  atmospheric  FEM  applications  is  still  in  its 
infancy  it  is  not  known  if  this  phenomenon  is  an  artifact  of  the  time 
stepping  scheme  employed,  the  boundary  conditions  imposed,  or  the  basic 
technique. 

A  simplistic,  second  order  diffusive  damping  term  was  included  in 
the  equations  of  motion.   To  keep  the  Galerkin  formulation  of  the 
diffusion  consistent  with  first  order  terms,  an  area  integration  was 
performed  and  Gauss'  divergence  theorem  was  applied.   Consider  the  diffu- 
sion term  in  the  zonal  equation  (III-2)  with  Galerkin' s  technique  applied: 


K    yV2uV.dA  =  K  f   V'VuV.dA 

ha     '     ha     : 


=  K   {yv.  Vu*n  dr  -J   W .  *VudA}  (IV-1). 


*■-■-      3  ~  3 


where  dr   is  differential  distance  along  the  path  of  integration  on  the 

outer  boundary,  K  is  the  horizontal  viscosity  coefficient  and  n   is 

H 

the  unit  normal  to  the  area  or  domain,  A. 

The  contour  integral  in  (IV-1)  vanishes  due  to  cyclic  continuity  and 
an  implied  -r—  =  0  .   The  diffusion  term  in  the  meridional  equation  is 
advanced  similarly. 
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V.   EOUATION  FORMATION 

As  an  example  of  the  use  of  the  Galerkin  technique  consider  equation 
(III-l)  .   A  consistent  form  for  (j>  was  assumed  and  substituted  into 
(II1-1) .   Then  an  area  integration  was  performed. 

Defining  the  inner  product  of  a  basis  function  g(x,y)   and  a  test 
function,  V(x,y)  as: 

/  /  gV  dxdy  =  <g,V> 
y  x 

An  application  of  the  FEM  to  the  continuity  equation  gives: 


<Hv>--<^  (u*>'  y>-<-h  M)- v> 


A  partial  integration  of  the  x-derivative  term  yields: 


8x 


■  (u<J>),  V>  =  /  -r-  (u(f)V)dA  -  /  u<J>  -*-  dA 


A 

But  the  first  term  on  the  right  hand  side  vanishes  because  of  cyclic 
continuity  in  x.   Similarly,  since  v  =  0  on  the  boundary,  the  correspond- 
ing term  in  the  partial  integration  of  the  y-derivative  term  also 
vanishes o 

Applying  Galerkin' s  linear  technique  it  is  assumed  $,  u,  and  v  are 
of  the  form 

<j>  =  Y-V. 
u  =  a.v. 

3    3 

v  =  B.v. 

3  3 
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which,  with  a  concurrent  test  function  assumption,  transforms  the  con- 
tinuity equation  to: 


Y  .  <V.  V.>  =  Y.  [<a,  V,  V.  V.  >  +  <R  V,  V.  V.  .>]  . 

'  3        3/  i     j    k  k  ],  ix      k  k  j,  i] 


where  the  dot  denotes  a  time  derivative. 

Using  a  centered  finite  time  difference  the  final  form  of  the 
equation  is: 


(Y.   -y.   )<v.  v.>  =  2At  Yn  [<a,nv,  v.  v.  >  +  <3,nv,v.  V.  >]   (V-l) 
3  3  3r     i         D    k  k  3t    ix      k  k  ]f  xx 


with  n  being  a  time  level. 

Equation  (V-l)  is  a  matrix  equation  of  the  form 


[A]  {x}   -   {b} 


where: 


[A]   =   <V.  V.> 
3'    i 


r  i      .   n+1     n-1, 
(x)   =  (y     .       -  Yj   ') 


and 


{b}    the  right-hand  side  of  (V-l) . 
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The  u  and  v  momentum  equations  were  similarly  advanced  including 

diffusion  and  become,  respectively: 

-a .n [<a,  nv  v .  ,v.>+<3,  nv,  v.  ,v.>] 
3    k  k  jx  1    k  k  jy      1 

(0tjn+1  '  ajn~1)<Vj'V  =  2At  (   +2fi  <sineo6knVk'Vi>"<YknVkx'Vi>1 

"Ku  ^'V  'V-  >+«*,  "v,  ,  V.  >] 
H    k  kx  ix    k  ky   ly 

0.n[<a1nv,v.  ,v.>+<ft  nv,  v.  ,vj 

3    k  k  jx  x    k  k  jy   i 


n+1     n-1.  __  _  ^    ..   /     __   .  0   n 


(6.     -  $.    )<V.,V.>  =  2At 
J       3      31  * 


+  2SJ<sineoak'vk/v.>+<YkVky,v.> 
■H    [<Sk\x'Vix>  +  <6k\y'Viy>: 


At  this  point  it  is  noted  that  all  the  inner  products  are  functions 
only  of  space  and  need  be  computed  but  once.   Also  useful  is  the  fact 
that  the  coefficient  matrix,  <V.,V.>  ,  is  identical  in  all  three  equa- 
tions.  Finally,  at  this  stage  all  equations  are  uncoupled  at  any  time 
step. 

As  is  apparent,  the  time  discretization  scheme  chosen  was  the  leap- 
frog or  centered  difference  method.   It  was  chosen  primarily  for  its 
simplicity.   Obvious  drawbacks  are  its  inherent  computational  mode  and 
the  somewhat  more  restrictive  maximum  time  step  than  some  semi-implicit 
schemes  may  allow.   The  latter  factor  was  expected  to  and  did  play  a 
major  role  in  the  amount  of  computational  time  involved  since  the 
smallest  grid  resolution  essentially  limits  the  maximum  allowable  time 
step.   It  was  expected,  based  on  the  work  of  Cullen  (1973)  using 
simple  linear  advection,  that  the  time  step  criterion  would  be  slightly 
more  restrictive  than  a  finite  difference  scheme  with  equivalent 
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resolution.   The  leapfrog  scheme  used  a  half -forward,  half  leapfrog  time 
step  to  start.   Also,  early  experiments  showed  strong  solution  separation, 
This  tendency  was  kept  in  check  with  a  periodic  restart,  identical  to 
the  initial  start,  being  applied  every  12  time  steps. 
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VI.   TREATMENT  OF  BOUNDARY  CONDITIONS 

The  cyclic  boundary  condition  is  naturally  included  in  the  system  by 
defining  how  the  nodal  points  are  connected.   Remaining  conditions  are 
those  imposed  at  the  northern  and  southern  walls.   The  geostrophic  condi- 
tion is  essentially  independent  of  time,  so,  to  fit  this  condition  into 
the  model,  boundary  equations  at  distinct  times,  n+1  and  n-1,  were  sub- 
tracted.  That  is: 

t|£  -2ftsin9  u)n+1  -  (|£  -2fisin6  u)n_1  =  0 
3y        o         dy        o 

The  Galerkin  formulation  gives: 


(Y.n+1-Y.n  1)  <v.  ,v.>  =  -f  [<a  n+1v.,v.>  -  <a  n~1v.,v.>]       (Vi-1) 
3  3  DY1      °   D    3  *      :    3      1 

with      f  =  2fisin0 
o         o 


The  rigid  wall  boundary  formulation  merely  required  setting  the 
Galerkin  approximation  for  meridional  wind  to  zero  at  the  i  boundary 
points  initially,  and  requiring  no  time  change  in  these  values: 

(3n+l   _  ^n-lj  <v      m   Q 
J         3  3   1 

These  boundary  conditions  were  applied  at  all  points,  i,  on  northern 
and  southern  boundaries. 

It  should  be  realized  that  the  geostrophic  boundary  formulations 
imply  that  a  different  coefficient  matrix  is  formed  than  that  shown  in 
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equations  (V-l,  V-2,  V-2) .   The  only  terms  in  these  matrices  which  are 
different  are  those  in  the  boundary  equations. 
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VII.   COMPUTATIONAL  TECHNIQUES 

A  primary  operational  disadvantage  of  any  finite  element  technique  is 
the  requirement  of  large  amounts  of  computer  storage.   This  results  from 
the  system  of  linear  equations.   The  problem  is  not  as  bad  as  it  first 
seems  since  in  any  one  equation  of  an  N  x  N  system  a  majority  of  coeffi- 
cient matrix  elements  are  zero.   The  problem  then  reduces  to  inverting  a 
system  with  a  very  sparse  coefficient  matrix.   If  one  can  utilize  an 
equation  solver  which  takes  advantage  of  this  property,  core  requirements 
can  be  minimized. 

Specifically,  for  the  ith  nodal  point,  the  number  of  non-zero  entries 
in  the  coefficient  matrix  for  equation  i  is  one  plus  the  number  of  elements 
surrounding  that  nodal  point.   Some  storage  schemes  take  advantage  of  the 
banded  nature  of  the  coefficient  matrix.   The  band  width  is  determined 
by  the  largest  different  between  the  number,  i,  of  the  ith  equation  and 
the  numbers  of  the  nodal  points  locally  connected  to  nodal  point  i.   Mini- 
mal bandwidth  is  thereby  gained  by  choosing  an  efficient  nodal  point 
numbering  scheme. 

It  was  decided  not  to  minimize  bandwidth  structure  for  any  cases  under 
consideration  since  the  time  involved  in  selecting  an  efficient  number 
scheme  for  differing  domains  possessing  cyclic  continuity  could  be  better 
spent.   A  disadvantage  of  a  banded  matrix  scheme  is  the  storage  of  some 
zero  values;  a  possible  problem  when  large  systems  are  to  be  solved.   This 
author's  unbanded  approach  has  the  disadvantage  of  forcing  the  use  of  an 
iterative  solution  technique  for  the  linear  system. 
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To  fully  compact  the  coefficient  matrix  a  novel  technique  was  used 
which  was  developed  by  Salinas  (private  communication) .   Essentially, 
only  non-zero  coefficient  matrix  entries  are  retained.   The  coefficient 
matrix  is  stored  as  a  linear  vector.   Three  vectors  herein  called  NUM, 
ISTART,  and  NAME,  are  required  to  use  the  method.   NUM  contains  for 
each  j  nodal  point  the  number  of  points  interacting  with  j  including 
itself.   NAME  gives  the  information  of  which  nodal  points  interact  with 
j  including  j.   ISTART  tells  where,  in  NAME,  the  numbers  of  the  points 
interacting  with  j  begin.   NUM  and  ISTART  are  both  dimensioned  N  for  a 
system  with  N  unknowns.   NAME  is  dimensioned  to  the  length  of  non-zero 
entries  in  the  coefficient  matrix.   That  is,  NAME'S  dimension  is  the 
sum  of  the  number  of  nodal  point  interactions. 

The  matrix  assemblage  required  is  mathematically  equivalent  to 
multiplying  each  basis  function  on  a  nodal  point  by  nodal  point  basis 
by  the  correct  test  functions,  and  then  integrating  over  the  portion  of 
the  domain  spanned  by  these  functions.   Graphically,  as  shown  in  Figure 
4,  the  basis  function  at  j  is  shown  outlined  in  heavy  black  over  a 
group  of  elements  around  nodal  point  j .   Portions  of  the  test  functions 
which  interact  with  this  basis  function  are  also  shown. 

Since,  in  the  case  of  linear  elements,  the  domain  integration  is 
analytically  obtained  over  an  element,  the  matrix  assemblage  of  an  inner 
product,  such  as  <V.,V.>  ,  is  done  most  efficiently  on  an  elementary 
level.   To  accomplish  this,  a  local  numbering  scheme  is  required  for  each 
element.   This  table  was  stored  in  the  two  dimensional  array,  NPTS  (N,3). 
NPTS  contains,  for  each  element,  e,  a  set  of  the  three  global  nodal 
points  which  define  e.   For  each  element  these  numbers  are  stored 
sequentially  according  to  number  of  the  second  dimension  of  NPTS.   The 
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Figure  4.   A  basis  function  and  portions  of  two  test  functions, 
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choice  of  the  point,  k,  to  start  numbering  e  was  arbitrary,  but,  once 

chosen,  the  remaining  two  points  were  numbered  in  sequence  proceeding 

in  a  positive  sense  around  e.   For  example  in  Figure  4  element  e  is 

labeled  NPTS  (e,l)  =  k,  NPTS  (e,2)  =  1,  and  NPTS  (e,3)  =  j. 

To  illustrate  matrix  assembly  using  an  element  by  element  technique 

consider  again  Figure  4.   The  computational  technique  required  is  for 

each  point  describing  e;  namely  k,  1,  and  j  ;  the  inner  product  value  be 

distributed  to  its  proper  location  in  the  coefficient  matrix,  <V  ,V.>  . 

Two  of  the  three  test  functions  over  e  are  drawn  labeled  t,  and  t,  . 

1      k 

Into  the  equation  for  point  k  a  contribution  from  the  interaction  between 
b   (not  shown)  and  t  must  be  added;  that  is,  <b  , t  >  (known  analytically) 

K  K.  K   K 

must  be  added  into  the  kth  element  of  equation  k.   Also  t  '  s  interaction 

with  b  must  be  entered  into  the  1th  element  of  equation  k.   Finally, 

t.'s  interaction  with  b,  must  be  added  into  the  jth  element  of  equation  k. 
j  k 

The  same  local  dispensing  of  interactions  must  be  carried  out  for  the 
remaining  nodal  points,  1  and  j  of  e.   This  procedure  when  done  for  all 
elements  in  the  domain  gives  <V.,V.>  or  the  coefficient  matrix  of  the 
equation. 

The  inner  products  used  may  either  be  computed  as  needed  during  the 
equations '  time  integration  or  computed  prior  to  time  integration  and 
stored  for  all  time.   The  latter  method  was  followed  since  the  original 
computer  system  used  was  in  IBM  360/67  which  is  relatively  slow  by  today's 
standards  but  was  able  to  handle  the  core  requirements  of  the  problems 
tested.  Also,  since  varying  size  elements  would  be  used  and  the  equations 
require  numerous  time  steps  for  even  a  six  hour  forecast,  an  inner  product 
storage  is  most  efficient. 
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VIII.   DOMAINS  CONSIDERED 

A.   BASIC  MESHES 

To  gain  experience  with  the  FEM,  numerous  grid  meshing  schemes  were 
considered.   The  first  problem  chosen  was  one  which  required  a  minimal 
amount  of  central  memory  storage  and  execution  time.   This  mesh  is  pre- 
sented in  Figure  5(a).   Actual  tests  were  conducted  on  grids  similar  to 
the  remaining  portions  of  Figure  5.   In  all  cases  the  domain  is  a  bounded 
channel,  cyclic  in  the  east-west  direction  with  a  width  and  height  of 
3503  km. 

Higher  resolution,  regular  grids  were  also  used.   One  of  these  is 
shown  in  Figure  5(b).   The  domain  was  divided  into  12  latitudinal  and 
longitudinal  subdivisions  (12x12) .   A  finer  resolution  used  on  this  domain 
had  twice  the  resolution  as  Figure  5(b),  (24x24),  and  is  not  shown.   A 
coarser  scheme  (not  shown)  had  six  latitudinal  and  longitudinal  sub- 
divisions, (6x6) . 

Irregular  grids  were  also  used.   Figure  5(c)  shows  the  first  of  these 
tested  wherein  a  central  region  in  the  x-direction  had  half  the  external 
mesh  size  covering  half  the  channel  length.   A  similar  grid  was  used  with 
a  varying  resolution  in  y.   These  two  were  combined  to  give  the  grid 
shown  in  Figure  5 (d) . 

Other  irregular  grids  were  used  in  which,  rather  than  abruptly  chang- 
ing the  resolution,  a  smooth  gradation  was  employed.   The  x  varying  ver- 
sion is  shown  in  Figure  5(e)  along  with  the  version  which  varies  in  x  and 
y,  Figure  5(f).   This  zonal  and  meridional  variation  was  made  with  no 
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particular  plan  other  than  a  smooth  transition  to  an  inner  region  of 
half  the  external  or  near  boundary  element  sizes. 

All  the  grids  so  far  mentioned  were  coded  to  allow  a  reversal  in 
diagonal  slant.   Intuitively,  from  examination  of  Figure  5,  it  would 
be  surprising  if  some  bias  was  not  introduced  when  domains  are  divided 
with  a  built  in  diagonal  preference.   These  same  grids  were  also  tested 
with  an  alternating  diagonal  slant  similar  to  that  shown  in  Figure  6. 

B.   FINAL  MESH 

In  order  to  study  a  relatively  small-scale  vortex,  a  meshing  scheme 
was  developed  whose  spatial  resolution  increased  radially  toward  the 
center  of  the  vortex.   The  grid  developed  is  shown  in  Figure  6.   The 
routine  used  to  locate  the  nodal  points  was  executed  independently  of 
the  time  integration  and  its  results  were  read  in  prior  to  overall  inte- 
gration.  Input  data  included  nodal  point  locations  and  the  local 
elementary  numbering  scheme,  NPTS. 

This  program  was  made  flexible  in  order  to  allow  for  a  varying  reso- 
lution over  an  even  number  of  latitude/ longitude  blocks  in  the  total 
domain.   Also  included  was  a  scheme  by  which  the  radial  resolution  on 
the  subdomain  could  be  varied  in  a  regular  manner.   This  procedure  involved 
decreasing  the  distance  between  nodal  points  along  the  diagonals  which 
extend  from  the  corner  points  of  the  subdomain  by  increasing  even  powers 
of  an  arbitrary  factor,  K.   This  is,  for  the  five  diagonal  divisions  shown 
of  length  d 


3 

d.  = 
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with  K  <  1  d  =  ungraded  domain  diagonal  length 

In  the  ungraded  region  the  diagonal  slants  were  varied  as  shown  to  minimize 

any  bias. 
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Figure   6.      Final  mesh. 
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IX.   INITIAL  CONDITIONS 

A.   SINUSOIDAL 

For  the  basic  testing  on  the  rectangularly  subdivided  grids  such  as 
shown  in  Figure  5,  simple  initial  conditions  were  used  which  consisted  of 
a  single  wave  in  the  x  direction  and  one  half  wave  in  the  y  direction. 
Such  initial  conditions  which  satisfy  the  non-linear  balance  equation  with 
f  constant  are  given  below: 


d)  =  f  A  sina,  cosa,  -  f  U(y-y  )  +  d) 
o      1     2    o     m     o 

—  ATT 

u  =  u  —  cosa,  cosa^  (IX-1) 

w      1     2 

2TTA   . 
v  =  -  — —  sina  sma 
X_        I      2. 


with         A  =  arbitrary  perturbation  amplitude 

W  =  channel  width  =  3500  km 

X  =  channel  length  =  3500  km 
L 

U  =  mean  zonal  wind  =  10  m  sec 

y   =  mid  latitude  value  of  y  =  W/2 
m 

(J)   =  mean  free  surface  goepotential  height  chosen  as 

4  2    -2 
4.9  x  10  m  sec 

Ot   =  TTy/W 

a0    =  27TX/X 
z  L 

Two  values  of  A  were  tested;  one  which  gave  a  maximum  perturbation  zonal 
wind  of  5.5  msec   and  another  for  0.55  msec   .   A  mean  zonal  flow  of 
10  msec   was  imposed  in  some  experiments  with  this  initial  condition. 
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The  initial  conditions  with  and  without  zonal  flow  are  presented  in 
Figures  7  and  8  respectively.   The  perturbation  amplitude  in  both  cases 
is  for  a  perturbed  zonal  wind  of  5.5  msec 


B.   VORTEX 

Analytic  initial  conditions  were  used  for  a  small  scale  perturbation. 
A  tangential  wind,  Vn,  was  defined  in  cylindrical  coordinates  as: 

r     2 

V  =  A  —  sin  (cr)         r  <  r  (IX-2) 

T     rQ  -  0 

V  =  0  r  >  r 
T  o 

with     r  =  radial  distance  from  center  of  perturbation 
r  =  maximum  radial  extent  of  perturbation 
c  =  7T/r0 

A  =  arbitrarily  chosen  amplitude  in  msec 
An  analytic  value  of  <J)*(r),  the  disturbance  portion  of  geopotential, 
can  be  obtained  from  the  gradient  wind  equation: 

2 

^3L  +  f  v  =  I*'  (ix-3) 

r     o  T   9r 


Integrating  (IX-3)  radially  outward  from  the  center  of  the  vortex 
(r=0)  one  obtains: 

(j)'  (r)  =  — j  {"jT"  r  "  j-[sin(2cr)  -  —  sin(4cr)  ] 

i 

=rr    [cos(2cr)  -  —  cos(4cr)]}  +         (IX-4) 

8c 

f  A  r  r    r    .  ,  _   .  1      .  _   ,  -i 

—  i- —  sin(2cr)  -  — -  cos(2cr) j 

r^  4    4c  .2 

0  8c 

+   constant  of  integration 
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Figure  7.   Initial  sinusoidal  conditions  with  U  =  10  msec 
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Figure  8.   Initial  sinusoidal  conditions  with  U  =  0. 
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The  amplitude  of  the  disturbance,  A,  was  chosen  to  give  a  desired 
maximum  tangential  wind  component  of  25  meters  per  second,  a  realistic 
value  for  tropical  models  at  500  mb.   The  mean  500  mb  geopotential 
height,  (J)  ,  was  added  to  this  (f> '  disturbance  field.   Also  added  was  a 
mean  slope  field,  similar  to  that  in  (IX-1) ,  to  allow  for  a  mean  zonal 
flow. 

The  radius  of  the  disturbance,  r  ,  was  chosen  so  that  the  vortex 
spaneed  the  graded  portion  of  the  domain  shown  in  Figure  6  equal  to  six 
largest  grid  distances;  thus  r  was  approximately  875  kilometers.   The 
initial  conditions  are  presented  in  Figure  9  with  mean  flow  component 
shown  only  in  the  (J)1  field. 
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Figure  9.   Initial  vortex  conditions  with  U  =  10  msec 
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X.       RESULTS 

A.       DISCUSSION 

1.  Notation 

The  use  of  various  meshes  on  the  same  domain  enabled  a  detailed 
study  of  the  basic  properties  of  the  FEM.   Obviously  the  order  of  case 
study  proceeded  from  the  coarsest,  6x6  grid  mesh  to  an  18x18  mesh  with 
varying  resolution  in  both  directions  to  the  24x24  mesh  which  had  the 
best  overall  resolution.   To  simplify  the  task  of  referencing  all  cases 
considered  (see  Table  I). 

All  the  above  cases  were  integrated  through  48  hours.   A  maximum 
stable  time  step  was  roughly  found  using  Cullen's  (1973)  work  as  a  guide. 
Cullen  states  that  for  simple,  two  dimensional,  linear  advection,  a  stable 
FEM  time  step  is  1//3  times  that  of  the  corresponding  finite  difference 
(leapfrog)  scheme.   For  the  12x12  mesh  with  a  mean  flow  of  10  meters  per 
second,  Cullen's  criterion  gives  a  stable  time  step  of  8.8  minutes.   A 
nine  minute  time  step  was  found  unstable  on  the  12x12  while  an  eight 
minute  time  step  was  stable.   This  verifies  Cullen's  analysis  since  a 
corresponding  stable  finite  difference  time  step  is  slightly  in  excess 
of  15  minutes.   Time  truncation  errors  were  not  of  primary  concern  in 
this  study,  so  close  to  maximum  time  steps  were  always  used. 

2.  Harmonic  Analysis 

A  harmonic  analysis  routine  of  the  geopotential  fields  was  used 
at  specified  times  in  the  integration.   This  analysis  was  accomplished 
along  constant  latitude  circles  for  wave  numbers  through  three,  six  and 
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TABLE  I 

Experiment  Abbreviations 
(all  valid  at  48  hours) 


6x6 

[SD] 

12x12 

[SD] 

18x12     [SD] 

12x18     [SD] 

6x6 

[SID] 

12x12 

[SND] 

18x12    [S**D] 

12xl8G[SD] 

6x6 

[S2D] 

12x12 

[SID] 

18xl2G[SD] 

12x18     [VD] 

12x12 

[S2D] 

18x12    [VD] 

12xl8G[V**D] 

12x12 

[S*D] 

18xl2G[VD] 

12xl8G[VD] 

12x12 

[S*ND] 

18xl2G[V**D] 

12x12 

[VD] 

18x18 

[SD] 

24x24    [SD] 

18x18 

[SND] 

24x24    [VD] 

18xl8G[SD] 

18xl8G[SND] 

18x18 

[VD] 

18x18 

[VND] 

18xl8G[VD] 

where 


6x6 
12x12 

12x18 


number  of  rectangular  subdomains  in  X-Y  directions 


24x24 

G 
S 

S* 

S** 

V 

V** 

1 

2 

D 
ND 


smoothly  graded  (as  opposed  to  abruptly  varying) 

sinusoidal  initial  conditions  with  maximum  u  perturba- 
tion of  5.5  msec--'-  and  U=10  msec" 

same  as  S  but  one-fifth  of  the  perturbation  amplitude 

same  as  S  but  U=0 

vortex  initial  conditions  with  maximum  tangential  wind, 
VT  =25  msec-1  and  U=10 

same  as  V  but  U=0 

diagonal  slant  from  upper  left  to  lower  right 

diagonal  slant  from  lower  left  to  upper  right 
(no  number  implies  alternating  slant) 

diffusion  included 

no  diffusion 
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12  depending  on  the  mesh.   Longitudinally  averaged  wave  amplitude  infor- 
mation was  obtained  and  used  in  data  reduction.   These  amplitude  data 
were  used  to  measure  generated  noise  and  damping  effects.   Phase  informa- 
tion was  used  to  obtain  wave  speed  by  computing  phase  changes  at  two  hour 
increments . 

The  spectral  analysis  scheme  required  equally  spaced  data  along 
a  latitude  circle.   Since  some  grids  tested  had  unevenly  spaced  nodal 
points,  the  interpolation  property  of  FEM  was  used  in  these  cases.   If  a 
value  of  a  field  variable,  z,  is  known  at  an  element's  vertices  located 
at  points  P.  (x. ,y . )  (i=l,2, 3) ;  then  the  following  is  true: 

z   =  a  +  bx,  +  cy 
z2  =  a  +  bx2  +  cy2 
z3  =  a  +  bx3  +  cy3 

Solving  this  3x3  system  for  a,b,  and  c;  one  can  obtain  z  at  any 
point  P  interior  to  the  element. 

3.   Matrix  Inversion  Technique 

A  straightforward,  but  time  consuming,  Gauss-Seidel  technique  was 
used  for  linear  system  solution  for  vectors  of  the  three  variables  in 
each  time  step.   It  was  noted  that  10  to  15  passes  were  necessary  for 
system  convergence  with  no  explicit  boundary  conditions.   Eighteen  to  23 
passes  were  required  when  boundary  conditions  were  imposed.   The  condi- 
tions were  explicitly  imposed  in  the  continuity  and  meridional  momentum 
equations.   Also  of  interest  is  the  fact  that  convergence  was  only 
possible, when  applying  boundary  conditions,  if  the  direction  of  G-S 
iteration  was  reversed  every  pass. 
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Convergence  was  considered  to  be  achieved  when  the  following 
criterion  was  met: 


n    n-l| 

i    i    max  over  i     ^  _      -6 
<  1  x  10 


l  '  max  over  i 


with:      n  being  a  G-S  iteration 

i  being  a  nodal  point  number 


x.  being  an  element  of  the  solution  vector,  {x} 


A  more  efficient  scheme  such  as  sequential  over-relaxation  is 
desirable  to  save  time.   Actual  project  completion  required  the  use  of 
many  meshes  all  of  which  probably  possess  different  optimum  over-relaxa- 
tion factors,  so  G-S  was  used  for  expediency. 
4.   Diffusion  Coefficients 

A  working  maximum  stable  diffusive  coefficient,  IC  ,  which  is 
equal  to  a  comparable  finite  difference  value  was  defined  as: 


k   =  Ax  /4At 
hm 


Using  this  as  a  maximum,  various  values  of  K  were  tried  on  the 

12x12  grid.   Values  ranging  from  .0005  to  .05  of  K   were  tried  and  the 

spectral  and  visual  results  compared.   The  ideal  was  to  use  as  little 

possible  diffusion  but  enough  to  control  generation  of  noise  in  higher 

wave  numbers.   A  subjective  value  was  obtained  of  .005  IC   which  was 

used  in  all  tests  including  diffusion.   The  size  of  K,   was  limited  by 

hm 

the  minimum  grid  separation  of  each  mesh.   Values  used  are  shown  below: 
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TIME  STEP 
(minutes) 

MINIMUM 

GRID  SPACING 
(km) 

Kh 
(m  sec   ) 

16 

580 

4.44x10 
2.22x10 

8 

290 

4 

145 

l.llxlO5 

0.5 

8.26 

(see  text) 

GRID 

6x6 

12x12 

18x12  &  18xl2G 

12x18  &  12xl8G 

18x18  &  18xl8G 

24x24 

Final  graded 


5,   Data  Reduction  and  Presentation 

The  following  noise  or  roughness  measure  was  defined: 

N 


£  i 


i=2 

Roughness  =  

ai 

with    a  =  zonally  averaged  amplitude  from  spectral 
analysis 

i  =  wave  number  index 

N  =  maximum  Fourier  wave  number  (3,  6,  or  12) 

A  measure  of  damping  was  used  which  was  the  ratio  of  amplitude 
of  any  Fourier  component  geopotential  perturbation  to  its  amplitude  at 
t=0.   This  was  used  to  determine  the  detrimental  effects  of  diffusion 
on  the  solutions. 

An  analytic  wave  speed  was  derived  from  an  average  propagation 
speed  determined  by  phase  angle  changes.   The  quasi-geostrophic  value 
of  the  phase  speed,  c,  is  given  by 


2      2     2 
(k    +  k    -  A  ) 


c  =  U 


x     X 


„  2    ,2    \2s 
(k   +  k   +  A  ) 

x     y 
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with  k  =  wave  number  in  x  direction 

x 

k  =  wave  number  in  y  direction 

y      9 

x2  =  f  V4> 

o  To 


d>  =  mean  free  surface  height 

-1 
U  =  10  msec 


This  formula  gives  phase  speeds  of  9.857,  9.958,  9.981,  and 
9.989  meters  per  second  for  wave  numbers  1,  2,  3,  and  4  respectively. 

Graphical  plots  were  produced  for  all  tests.   Those  not  dis- 
cussed as  part  of  the  results  are  presented  in  Figures  35  through  54 
inclusive.   The  results  for  the  entire  length  of  the  cyclically  con- 
nected channel  were  not  plotted.   Each  plot  extends  to  within  one 
regular  grid  increment  of  the  boundary  where  the  mesh  is  connected. 
The  maps  of  <J)   and  u  have  the  domain  mean  removed.   The  nodal  points 
are  marked  with  small  crosses. 

B.   SINUSOIDAL  RESULTS 

The  effects  of  diagonal  slant  can  be  seen  in  Figures  10  and  11. 
The  most  noticeable  difference  exists  in  Figure  10  where  the  mirror 
image  effect  of  diagonal  bias  shows  itself  in  the  latitudinal  direction. 
The  cases  shown  possessing  an  alternating  slant  form  a  rough  average 
between  the  opposite  slant  cases.   Figure  11,  which  shows  the  $      ampli- 
tude, gives  a  slight  hint  of  bias  effects,  but  the  tendencies  are  less 
pronounced.   Wave  speed  data  for  similar  cases  on  the  24x24  grid  were 
not  obtained,  but  one  would  expect  less  bias  as  the  element  size 
approaches  zero.   A  plot  of  wave  number  one's  amplitude  versus  latitude 
for  the  12x12  mesh  showed  no  significant  variation  for  the  three  differ- 
ent diagonalization  methods,  so  apparently  the  effect  does  diminish 
with  finer  resolution. 
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Figure  10.   Latitudinal  variation  of  wave  speed  (48  hour  time  averaged) 
for  four  grids. 
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The  latitudinal  variation  of  phase  propagation  of  various  grids  is 
shown  in  Figure  12.   In  general  all  phase  speeds  are  greater  than 
analytic,  but  not  significantly  so.   It  should  be  noted  that  the  scale 
of  the  abscissa  in  Figure  12  covers  a  very  narrow  range  of  phase  speeds; 
even  the  coarsest  grid  is  within  four  per  cent  of  the  true  phase  speed. 
The  phase  velocities  depart  most  from  analytic  near  the  boundaries, 
where  the  amplitude  is  small.   A  somewhat  anomalous  result  is  the 
better  treatment  of  propagation  by  the  12x12  as  opposed  to  the  24x24 
mesh.   This  is  most  apparent  below  the  mid-channel  latitude.   No  explana- 
tion is  offered  for  this  condition. 

The  necessity  of  some  spatial  smoothing  is  apparent  when  one  compares 
Figures  13  and  14;  Figures  15  and  16;  Figures  17  and  18;  and  Figures  19 
and  20.   Even  if  i   were  not  for  this  unwanted  noise,  large  phase  errors 
result  without  smoothing,  as  Figure  21  suggests.   Even  by  decreasing  the 
initial  disturbance  amplitude  by  a  factor  of  five,  the  phase  errors  are 
almost  identical  with  unsmoothed  tests. 

The  fine,  course,  and  intermediate  resolution  grid  results  are  shown 
for  their  roughness  properties  in  Figure  22.   Even  with  four  times  the 
damping,  the  6x6  mesh  cannot  compete  with  the  24x24  mesh  in  this  compari- 
son.  The  6x6  results  become  even  less  desirable  when  one  sees  that  the 
wave  amplitude  is  cut  more  than  in  half  (Figure  23) .   An  enlarged  view 
of  Figure  22  is  presented  in  Figure  24.   This  diagram  should  show  if  any 
significant  improvement  can  be  gained  by  using  a  grid  which  is  inter- 
mediate between  the  12x12  and  24x24  mesh,  thus  increasing  accuracy  with 
less  computational  sacrifices.   It  should  be  noted  that  all  these  inter- 
mediate grids  possess  one  half  the  diffusion  of  the  12x12  mesh. 
Apparently  two  of  the  grids  tried  do  offer  the  desired  quality,  the  12x18 
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Figure  12.   Latitudinal  variation  of  wave  speed  (48  hour 
time  averaged) . 
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Figure   13.      12x12    [SD]    at   t=43   hours. 
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Figure   14.      12x12  [SND]    at   t=48  hours. 
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Figure   15.      12x12    [S*D]    at   t=48   hours 
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Figure   16.      12x12    [S*ND]    at   t=48   hours. 
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Figure  17.   18x18  [SD]  at  t=48  hours. 
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Figure   18.      18x18    [SND]    at   t=48  hours. 
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Figure  19.   18xl8G  [SD]  at  t=48  hours 
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Figure  20.   18xl8G  [SND]  at  t=48  hours. 
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Figure  21.   Latitudinal  variation  of  wave  speed  (48  hour 
time  averaged)  with  no  diffusion. 
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Figure  23.   Damping  characteristics  of  various  grids  at 
t=48  hours. 
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and  12xl8G  meshes.   Both  of  these  meshes  have  better  resolution  over 
the  region  where  the  disturbance  is  strongest.   It  was  hoped  that  the 
grid  with  a  graded  rather  than  abrupt  change  in  the  y  direction  would 
be  better.   However,  it  appears  the  graded  version  (12xl8G)  is  only 
better  at  latitudes  south  of  15  North.   The  18x12  and  18xl2G  meshes 
are  generally  worse  than  the  12x12  counterpart  although  the  18xl2G  is 
the  better  of  the  two.   This  perhaps  is  the  result  of  the  disturbance 
being  advected  across  this  graded  or  abruptly  changing  region  of  the 
two  grids.   One  would,  however,  hope  that  the  18xl2G  would  show  better 
noise  properties  than  the  12x12.   It  is  plausible  that  the  element  size 
change  is  still  too  abrupt  to  describe  the  scale  of  motion  used  with 
an  advecting  velocity.   The  remaining  two  cases  to  be  identified  are 
the  18x18  and  18xl8G.   The  18x18  case  shows  an  oscillating  effect  pro- 
ceeding northward.   This  could  be  explained  as  some  pseudo-boundary 
influence  induced  by  the  abrupt  element  size  change.   The  18xl8G  handles 
the  situation  somewhat  more  effectively.   The  oscillation  is  still 
present,  slightly  reduced,  but  the  roughness  is  significantly  less  than 
the  18x18  case  and  is  better  than  the  12x12  case  at  most  latitudes. 

The  spectral  distribution  of  wave  energy  at  48  hours  is  shown  in 

2   -2 
Figure  25.   Some  grids  had  disturbance  amplitudes  of  less  than  1  m  sec 

(gpm)  and  therefore  don't  appear  on  this  logarithmic  plot.   As  expected, 

the  24x24  has  the  least  tendency  for  energy  transferral  to  higher  wave 

numbers.   Of  interest  is  the  large  amount  of  energy  exchanged  by  the 

18x12  grid.   It  has  the  largest  energy  levels  in  all  higher  wave  numbers 

save  wave  number  four,  where  its  smoothly  graded  companion  (12xl8G)  is 

slightly  higher.   Relating  the  two  18x18  grids,  one  can  see  an  expected 

improvement  when  a  graded  version  of  the  same  grid  is  used;  a  useful 

feature  which  is  not  as  evident  at  latitudes  different  from  15. 
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Figure  25.   Spectral  characteristics  with  sinusoidal 
conditions  at  t=48  hours. 
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As  an  interesting  sidelight,  consider  Figure  26.   The  gravity  wave 
fluctuations  due  to  slight  imbalances  are  indicated  on  this  diagram. 
The  free  surface  gravity  waves  allowed  in  this  case  have  an  analytic 
period  of  4.4  hours  calculated  for  a  wave  length  equal  to  that  of  the 
channel.   This  agrees  well  with  the  short  period  fluctuations  shown  in 
wave  speed . 

C.   VORTEX  RESULTS 

The  gravity  wave  influence  on  geopotential  height  variations  is  also 
implied  in  Figure  27.   The  diffusion  is  also  evident  as  the  stationary 
vortex  fills  through  48  hours.   The  gravity  wave  amplitude  in  this  case 
is  larger  initially  than  with  the  sinusoidal  case  because  of  the  larger 
truncation  errors  in  representing  the  initial  state.   Rapid  filling  oc- 
curs in  the  first  18  hours  due  to  diffusion. 

Again  diffusion  is  necessary.   To  verify  this  compare  Figures  28 
and  29.   The  undesirable  requirement  of  significant  diffusion  at  500  mb 
is  obvious,  but  in  the  vortex  cases  the  diffusion  imposed  significantly 
damps  the  disturbance.   Further  experiments  should  verify  if  the  scheme 
could  survive  with  less  diffusion  of  the  type  used. 

Figure  30  shows  a  view  of  damping  for  grids  other  than  the  6x6  and 
24x24.   Eighty  eight  per  cent  of  the  original  mid-latitude  perturbation 
is  contained  in  waves  one  through  three,  so  these  should  be  of  primary 
concern.   All  grids  in  this  figure  except  the  12x12  have  identical  damp- 
ing coefficients  and  therefore  should  treat  up  to  wave  number  three 
similarly.   This  is  the  case  except  both  12x18  grids  don't  damp  wave 
number  three  as  much.   Since  this  Laplace  type  filter  is  not  very  wave 
number  selective,  one  would  expect  higher  diffusion  of  longer  waves. 
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wave  speed  fluctuations. 
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Figure  28.   18x18  [VD]  at  t=48  hours, 
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Figure   29,      18x18    [VND]    at   t=48   hours. 
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Figure  30.   Damping  characteristics  with  vortex  at  t=48  hours, 
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It  appears  that  both  12^:18  meshes  are  not  consistent  with  this  assump- 
tion.  Both  damp  wave  number  three  about  the  same  as  two  but  treat  wave 
numbers  four  and  five  in  opposite  senses.   The  12xl8G  grid  also  shows 
the  same  tendency  as  the  12xl8G  though  it  damps  wave  number  three  more. 
The  18x12  and  18xl2G  grids  also  show  similar  treatment  of  waves  four 
and  five.   The  large  damping  of  wave  five  relative  to  four  can  be  ex- 
plained in  part  by  the  fact  five's  initial  amplitude  is  about  two  and 
one  half  times  that  of  four.   The  large  damping  of  wave  two  by  the  12x12 
mesh  is  a  serious  defect  of  the  required  amount  of  coarse  grid  diffu- 
sion.  The  actual  spectral  results  of  these  grids  are  shown  in  Figure 
31. 

Figures  32  and  33  show  the  final  distributions  of  amplitude  for  waves 
one  and  two,  respectively.   The  large  damping  may  make  it  hard  to  dis- 
tinguish between  results  for  grids  with  greater  resolution  than  the 
12x12.   If  this  is  but  a  secondary  effect  then  one  could  say  that  extend- 
ing the  entire  domain  resolution  to  that  of  the  24x24  is  not  necessary. 

A  comparison  of  mid  channel  wave  speeds  shows  all  grids  studied  give 
very  accurate  results  for  wave  numbers  one  and  two.   Phase  propagation 
of  wave  number  three  is  handled  poorly  by  the  12x12  grid  and  is  270% 
in  error  for  the  12x12  case  and  70%  in  error  for  the  18x12  cases.   All 
grids  but  the  24x24  are  grossly  in  error  for  wave  number  four;  at  best 
175%  error  (18x18  grids)  and  at  worst  240%  error  for  the  18xl2G  case. 
The  24x24  case  is  25%  in  error  for  wave  number  four.   All  per  cent  errors 
mentioned  are  over  estimates  by  the  grids  in  wave  speed.   Phase  predic- 
tion accuracy  rapidly  deteriorates  for  higher  wave  numbers.   Propagation 
errors  for  large  time  integrations  are  contaminated  by  artificial  noise 
being  filtered  to  higher  wave  numbers. 
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Figure   31.      Spectral   characteristics  with  vortex   conditions   at 
t=48  hours. 
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Figure  33.   Latitudinal  variation  of  geopotential  for  wave  number 
two  at  t=48  hours. 
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Studies  with  the  highly  variable  grid  shown  in  Figure  6  were  cur- 
tailed.  It  was  therefore  not  possible  to  fully  investigate  the  movement 
of  the  grid  with  the  vortex.   A  few  experiments  were  carried  out  using 
both  the  vortex  and  single  harmonic  wave.   The  experiments  with  simpler 
initial  conditions  are  not  presented.   A  study  was  made  of  the  vortex 
with  no  mean  flow.   The  results  of  this  investigation  are  shown  in 
Figure  34. 

Basic  tests  were  done  with  a  constant  tC  over  the  domain.   The  value 
for  this  coefficient  was  again  selected  based  on  the  smallest  grid  in- 
crement.  On  this  grid,  this  value  was  very  small,  so  small  in  fact 
that  another  formulation  of  diffusion  was  required. 

A  more  involved  diffusion  formulation  was  tried  which  assumed  a 
diffusion  coefficient  constant  over  each  element.   A  vector  type  diffu- 
sion was  then  applied.   Simply  stated,  the  diffusion  over  an  element 
was  derived  using: 

A 
ne    *   A   nm 

with:      K,   =  diffusion  on  an  element,  e 
he 

A   =  area  of  element  e  or  element  1 

K,  =  diffusion  coefficient  based  on  minimum  grid 
spacing,  associated  with  elements  1  through 
24  (elements  surrounding  the  nodal  point  in 
the  center  of  diagram  Figure  6) . 

The  plots  in  Figure  34  show  that  an  initial  roughness  developed 
evident  at  six  hours  due  to  initial  gradient  imbalance  which  at  later 
times  was  minimized  as  smoothing  and  geostrophic  adjustment  effects 
dominateo 
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Figure  34.   Geopotential  field  for  stationary  vortex  on  final  mesh. 
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This  final  case  was  integrated  on  both  an  IBM  360/67  machine  and  a 
CYBER  17  5  machine.   The  CYBER  machine  is  twice  as  accurate  so  machine 

truncation  effects  were  minimized.   A  comparison  of  spectral  results 

« 
from  both  machines  at  t=12  hours  showed  little  difference,  so  machine 

I 
truncation  is  not  of  concern. 
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Figure  35.   6x6  [SD]  at  t=48  hours. 
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Figure  36.   6x6  [SID]  at  t=48  hours. 
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Figure   37.      6x6    [S2D]    at   t=48   hours. 
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Figure  38.   18x12  [S**D]  at  t=48  hours. 
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Figure   39.      12x12    [SID]    at   t=48   hours. 
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Figure  40.   12x12  [S2d]  at  t=48  hours. 
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Figure  41.   18x12  [SD]  at  t=48  hours. 
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Figure  42.   18xl2G  [SD]  at  t=48  hours. 
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Figure  43.   12x18  [SD]  at  t=48  hours. 
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F.igure   44.      12xl8G    [SD]    at   t=48  hours. 
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Figure  45.   24x24  [SD]  at  t=48  hours. 


3360 


89 


a 

■  — ^-^ — — — .  f 

m 

~*^-^v    ~~~      *  -^ 

o 

^^—~ 

,' 

^N 

"-v^              

~                       ^ 

vK~- — ~_: 

. 

*-r 

S 

J9 

y^ - ' ' 

« 

..--■     -    '    _.       - 

>^_ 

,  -^ 

-r^            ^-^i 

—   -■    — —  .. 

" ~^.                     ^"~--~- 

~^^~--. 

"^~"     -                       "^«k. 

. 

-- 

'         "            -.                                             '---. 

v       V  ■ 


•  .  ' 


■BSf. 


^ 


,/ 


— Lll ^^.. 


161X3 

x(icm) 


\* 


\ 


-^ 


^2CT 


ML' 


•i'<S 


V 


<P   Disturbance    (gpm) 


g 

in 


>)■ 


■ 


.) 


) )  )    . 


—-c   .■■ 


'  y 


/ 


1600 

xOcm) 


r" 


--') 


"3200 


u    Disturbance  (m/sec) 


Figure   46.      12x12    [VD]    at   t=48   hours, 
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figure  47.   18xl2G  [V**D]  at  t=48  hours 


91 


<P   Disturbance    (gpm) 


T396" 

xOcm) 


u   Disturbance  (m/sec) 


figure   48.      18x12    [VD]    at   t=48   hours 
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Figure  49.   18xl2G  [VD]  at  t=48  hours. 
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figure  50.   12x18  [V**D]  at  t=48  hours 
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Figure  51.   12x18  [VD]  at  t=48  hours. 
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Figure  52.   12xl8G  [VD]  at  t=48  hours. 
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Figure  53.   18xl8G  [VD]  at  t=48  hours, 
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Figure   54.      24x24    [VD]    at   t=48  hours, 
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XI.   CONCLUSIONS 

In  this  research  a  finite  element  model  has  been  developed  for  the 
barotropic  primitive  equations.   The  motion  is  confined  in  a  periodic 
channel  between  rigid  walls.   The  model  was  tested  with  two  sets  of  analy- 
tic initial  conditions. 

The  experiments  with  a  sinusoidal  wave  showed  excellent  phase  pro- 
pagation when  compared  to  the  quasi-geostrophic  analytic  solution  and 
when  some  diffusion  was  present.   Even  the  6x6  grid  gave  an  error  of  only 
about  4  per  cent  which  is  no  larger  than  the  error  obtained  with  a  fourth 
order  finite  difference  model.   All  grids  produced  computational  noise 
without  diffusion  which  made  the  forecasts  unrealistic  at  t=48  hours. 

In  order  to  control  this  noise  Laplacian  diffusion  terms  were  added  to 

i 
the  momentum  equations.   The  low  resolution  experiments  required  a  much 

higher  diffusion  coefficient  than  the  high  resolution  experiments.   This 

lead  to  heavy  damping  of  the  original  wave  for  low  resolution  and  little 

damping  for  high  resolution. 

The  same  initial  disturbance  was  tested  with  grids  which  contained 
an  abrupt  change  in  element  size  and  grids  which  contained  a  gradual 
variation  in  element  size.   These  variations  in  element  size  did  generate 
some  noise,  but  this  noise  was  handled  by  the  previously  introduced 
diffusion.   The  smoothly  varying  elements  gave  somewhat  better  solutions. 

The  model  was  also  tested  with  an  isolated  vortex  with  the  initial 
pressure  determined  analytically  from  the  gradient  wind  equation.   The 
experiments  with  a  uniform  element  size  displayed  gravity  waves  which 
were  excited  by  the  finite  element  imbalance  in  the  initial  conditions. 
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The  required  diffusion  damped  and  enlarged  the  initial  vortex.   The  vortex 
experiments  with  a  variable  element  size  showed  much  better  results.   This 
was  a  result  of  the  higher  resolution  around  the  vortex  and  the  lower 
diffusion.   Some  experiments  were  performed  with  a  continuously  varying 
element  formulation  which  had  very  small  central  elements.   This  procedure 
was  very  uneconomical  because  a  very  small  time  step  had  to  be  used 
throughout  the  domain  (see  Table  II) . 

This  study  has  verified  the  previously  discussed  phase  accuracy  of 
the  finite  element  method.  However,  the  noise  generation  was  larger  than 
expected.   It  is  not  known  whether  or  not  some  of  the  noise  arises  from 
imperfect  boundary  conditions  at  the  walls.   Cullen  (1974,  1976)  has 
discussed  the  necessity  of  high  order  smoothing  with  finite  element  models, 
In  any  case  the  model  developed  in  this  study  would  be  benefited  by  a 
high-order  spatial  filter  which  would  control  the  noise  while  not  affect- 
ing the  larger  scale  features. 

The  main  problem  with  the  use  of  variable  element  size  is  the  require- 
ment that  the  same  time  step  be  used  everywhere.   Finite  difference 
models  need  use  small  time  step  only  in  the  high  resolution  regions. 
Until  this  problem  is  solved  it  would  appear  that  the  element  size  should 
not  be  varied  by  a  large  amount.   But  since  the  FEM  is  more  accurate 
perhaps  tropical  storm  forecasts  can  be  made  without  excessively  small 
element  sizes.   The  model  developed  in  this  paper  requires  testing  with 
a  moving  element  system.   Further  work  is  also  required  to  make  the  model 
more  efficient. 
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TABLE  II 

48  Hour  PROGNOSIS 

IBM  360/67  TIME  AND  STORAGE  REQUIREMENTS 

(Fortran  H  compiler) 


AVERAGE 

GRID 

EXECUTION  TIME 

6x6 

3  minutes 

12x12 

22  minutes 

12x18 

1  hour  20  minutes 

18x12 

18x18 

1  hour  50  minutes 

24x24 

3  hours  40  minutes 

Final 

10  hours  (IBM) 

Graded 

1  hour  20  minutes 

Grid 

STORAGE  REQUIREMENTS 
(IN  BYTES,  1  BYTE  =  4  WORDS) 

49K 
129K 

188K 

•   261K 
44  2K 


168K  (IBM) 
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